function DescriptivesPaper(opt)

[~,Action,~,AgeState,~,~,DS] = prepdataM10(opt);


%% Display table of descriptive statistics (Article Table 1)

% Mapbiomas land use classification:
LU = nan(size(DS,1),1);
LU(DS.mb_uso_detalhe == 3 | DS.mb_uso_detalhe == 9)  = 1; % Forests
LU(DS.mb_uso_detalhe == 4 | DS.mb_uso_detalhe == 12) = 2; % Savanna (including grasslands)
LU(DS.mb_uso_detalhe == 15)  = 3; % Pasture
LU(DS.mb_uso_detalhe >= 19 & DS.mb_uso_detalhe <= 21) = 4; % Agriculture
DS.pasture_dummy = (LU == 3);
DS.crop_dummy = (LU == 4);
DS.forest_dummy = (LU == 1);
DS.savanna_dummy = (LU == 2);

disp(' ')
disp(' Descriptive stats for fields characteristics ')

% First do Sugarcane road distance (because it expans several columns):
SRDloc = [16:25];
SRD = double(DS(:,SRDloc));
SRD = SRD(SRD>0)/1000;
fprintf(1,'%20s ',   'SRD');
fprintf(1,'%12.3f ', nanmean(SRD));
fprintf(1,'%12.3f ' , nanstd(SRD) );
fprintf(1,'%12.3f ' , min(SRD) );
fprintf(1,'%12.3f ' , max(SRD) );
fprintf(1,'%12.3f ' , quantile(SRD,.2));
fprintf(1,'%12.3f ' , quantile(SRD,.5));
fprintf(1,'%12.3f ' , quantile(SRD,.8));
fprintf(1,'\n');
clear SRD
FieldsCharacteristics = [58:60, 73, 38, 53, 74:77];

for ch = FieldsCharacteristics
    fprintf(1,'%20s ',   char(DS.Properties.VarNames(ch)));
    fprintf(1,'%12.3f ' , nanmean(double(DS(:,ch))) );
    fprintf(1,'%12.3f ' , nanstd(double(DS(:,ch))) );
    fprintf(1,'%12.3f ' , min(double(DS(:,ch))) );
    fprintf(1,'%12.3f ' , max(double(DS(:,ch))));
    fprintf(1,'%12.3f ' , quantile(double(DS(:,ch)),.2));
    fprintf(1,'%12.3f ' , quantile(double(DS(:,ch)),.5));
    fprintf(1,'%12.3f ' , quantile(double(DS(:,ch)),.8));
    fprintf(1,'\n');
end


%% Conditional probability of sugarcane expansion (Article Figure 3)
fig = 1;
Nboot = 200; % Number of bootstrap repetitions for CI in figure

disp(' Probability of expansion v. distance from existing sugarcane fields ')
fig = fig + 1;
figure(fig)

SugRoadDist = double(DS(:,16:25)); Xvar = []; XvarLagged = []; XvarLagged2 = []; Yvar = [];
for t=7:10
    NotSugBeg = (AgeState(:,t) == 1);
    SugarExpansionYear = (AgeState(NotSugBeg,t) == 1 & Action(NotSugBeg,t) == 1);
    SugRoadDistYear = SugRoadDist(NotSugBeg,t-1);
    SugRoadDistYearLagged = SugRoadDist(NotSugBeg,t-3);
    SugRoadDistYearLagged2 = SugRoadDist(NotSugBeg,t-5);
    Xvar = [Xvar;SugRoadDistYear];
    Yvar = [Yvar;SugarExpansionYear];
end
Nq = 100; LowExtreme = .01; HighExtreme = .36;
Q  = linspace(quantile(Xvar,LowExtreme),quantile(Xvar,HighExtreme),Nq);
[Yhat] = KernelRegression(Yvar,Xvar,Q');
n = length(Yvar);
disp(' Bootstraping for CI:')
for b=1:Nboot
    sel     = ceil(n*unifrnd(0,1,n,1)); % bootrap randomization
    Yb = Yvar(sel); Xb = Xvar(sel); % Select bootstrap data
    [Ypredb(:,b)] = KernelRegression(Yb,Xb,Q');
    ProgressBar(b,Nboot)
end
Yquant = quantile(Ypredb',[0.025; 0.975]);
plot(Q/1000,Yhat,'LineWidth',1.5,'LineStyle','-','Color','black')
hold on
plot(Q/1000,Yquant,'LineWidth',1,'LineStyle',':','Color','black')
ylabel('Probability of sugarcane planting')
xlabel('Effective road distance to closest sugarcane field (km on highway)')
clear Yb
saveas(gcf, 'Figures/sugar_distance', 'eps')

fig = fig + 1;
figure(fig)

disp(' Probability of expansion v. Sugarcane yield ')
NotSugBeg = (AgeState(:,3) == 1);
n        = sum(NotSugBeg);
Xvar     = DS.gaez_aecol_suc(NotSugBeg);
SugarExpansion = any(AgeState(NotSugBeg,:) == 1 & Action(NotSugBeg,:) == 1,2);
Nq = 100; LowExtreme = .03; HighExtreme = .97;
Q  = linspace(quantile(Xvar,LowExtreme),quantile(Xvar,HighExtreme),Nq);
[Yhat] = KernelRegression(SugarExpansion,Xvar,Q');
disp('Bootstraping for CI:')
% Bootstrap for CI:
for b=1:Nboot
    sel     = ceil(n*unifrnd(0,1,n,1)); % bootrap randomization
    Sb = SugarExpansion(sel); Xb = Xvar(sel); % Select bootstrap data
    [Yb(:,b)] = KernelRegression(Sb,Xb,Q');
    ProgressBar(b,Nboot)
end
Yquant = quantile(Yb',[0.025; 0.975]);
plot(Q,Yhat','LineWidth',1.5,'LineStyle','-','Color','black')
hold on
plot(Q,Yquant,'LineWidth',1,'LineStyle',':','Color','black')
ylabel('Probability of sugarcane planting')
xlabel('Sugarcane yield (DW kg/ha)')
saveas(gcf, 'Figures/sugar_potyield', 'eps')

%% Replanting age by cohort (Article Figure 2)

% Replanting probabilities by cohort:
figure(fig); fig = fig + 1;
cohorts_used = 2:6; 
c_idx = 1;
ReP = nan(11,length(cohorts_used));
for c = cohorts_used
    IsCohort  = (AgeState(:,c)==2);
    %CohortAge = AgeState(IsCohort,c:end);
    CohortAct = Action(IsCohort,c:end);
    Replant = (CohortAct == 1);
    age = 1;
    for t = c:10
        ReP(age,c_idx) =  mean(Replant(:,age));
        Replant(CohortAct(:,age)==1 | CohortAct(:,age)== -1,age+1:end) = 0; % If replanted at that age, treat all future values as 0;
        age = age + 1;
    end
    c_idx = c_idx + 1;
end
%plot(ReP/0.7528,'LineWidth',2)
hold on
plot(ReP(:,1),'k-','LineWidth',2)
plot(ReP(:,2),'k--','LineWidth',2)
plot(ReP(:,3),'k-.','LineWidth',2)
plot(ReP(:,4),'k:','LineWidth',2)
plot(ReP(:,5),'k-o','LineWidth',2)
legend('2003', '2004', '2005', '2006', '2007', '2008', '2009')
xlabel('Age')
ylabel('Fraction replanting')
hold off

% Cumulative replanting probabilities by cohort:
figure(fig); fig = fig + 1;
RePCDF = nan(size(ReP));
for a = 1:size(ReP,1)
    RePCDF(a,:) = sum(ReP(1:a,:),1);
end
hold on
plot(RePCDF(:,1),'k-','LineWidth',2)
plot(RePCDF(:,2),'k--','LineWidth',2)
plot(RePCDF(:,3),'k-.','LineWidth',2)
plot(RePCDF(:,4),'k:','LineWidth',2)
plot(RePCDF(:,5),'k-o','LineWidth',2)
legend('2003', '2004', '2005', '2006', '2007', '2008', '2009')
xlabel('Age')
ylabel('Fraction replanting')
hold off

end
